############################
#
# Addictive Behaviors Archive
# B. W. Campbell
# Last Updated:  01/15/2018
#
############################
# Top Matter
############

rm(list = ls())
options(stringsAsFactors = FALSE)

source("GraduationDataPrep.R")

library(tnam)

##############
# Affirmations
##############

Affirmations_1 <- tnam(success_list ~
                         #netlag(success_list, correct_adj, pathdist = 1) +
                         netlag(success_list, affirm_adj, pathdist = 1),
                       #netlag(success_list, correct_adj, pathdist = 2) +
                       #netlag(success_list, affirm_adj, pathdist = 2) +
                       #weightlag(success_list, correct_adj) +
                       #   weightlag(success_list, affirm_adj) +
                       # centrality(correct_adj, type = c("indegree")) +
                       #  centrality(affirm_adj, type = c("indegree")) +
                       #  covariate(reciprocity_rate_list)+
                       #covariate(days_list),
                       #  covariate(age_list) +
                       #  covariate(black_list) +
                       #  covariate(lsiExit_list),
                       re.node = FALSE,
                       time.linear = FALSE,
                       family = binomial(link="logit"))

Affirmations_2 <- tnam(success_list ~
                         #netlag(success_list, correct_adj, pathdist = 1) +
                         netlag(success_list, affirm_adj, pathdist = 1) +
                         #netlag(success_list, correct_adj, pathdist = 2) +
                         netlag(success_list, affirm_adj, pathdist = 2),
                       #weightlag(success_list, correct_adj) +
                       #   weightlag(success_list, affirm_adj) +
                       # centrality(correct_adj, type = c("indegree")) +
                       #  centrality(affirm_adj, type = c("indegree")) +
                       #  covariate(reciprocity_rate_list)+
                       #   covariate(days_list)
                       #  covariate(age_list) +
                       #  covariate(black_list) +
                       #  covariate(lsiExit_list),
                       re.node = FALSE,
                       time.linear = FALSE,
                       family = binomial(link="logit"))

Affirmations_3 <- tnam(success_list ~
                         #netlag(success_list, correct_adj, pathdist = 1) +
                         netlag(success_list, affirm_adj, pathdist = 1) +
                         #netlag(success_list, correct_adj, pathdist = 2) +
                         netlag(success_list, affirm_adj, pathdist = 2) +
                         #weightlag(success_list, correct_adj) +
                         #   weightlag(success_list, affirm_adj) +
                         # centrality(correct_adj, type = c("indegree")) +
                         centrality(affirm_adj, type = c("indegree")),
                       #  covariate(reciprocity_rate_list)+
                       #   covariate(days_list)
                       #  covariate(age_list) +
                       #  covariate(black_list) +
                       #  covariate(lsiExit_list),
                       re.node = FALSE,
                       time.linear = FALSE,
                       family = binomial(link="logit"))

Affirmations_4 <- tnam(success_list ~
                         #netlag(success_list, correct_adj, pathdist = 1) +
                         netlag(success_list, affirm_adj, pathdist = 1) +
                         #netlag(success_list, correct_adj, pathdist = 2) +
                         netlag(success_list, affirm_adj, pathdist = 2) +
                         #weightlag(success_list, correct_adj) +
                         #   weightlag(success_list, affirm_adj) +
                         # centrality(correct_adj, type = c("indegree")) +
                         centrality(affirm_adj, type = c("indegree")) +
                         covariate(reciprocity_rate_list, coefname = "Reciprocity Rate"),
                       #   covariate(days_list)
                       #  covariate(age_list) +
                       #  covariate(black_list) +
                       #  covariate(lsiExit_list),
                       re.node = FALSE,
                       time.linear = FALSE,
                       family = binomial(link="logit"))

Affirmations_5 <- tnam(success_list ~
                         #netlag(success_list, correct_adj, pathdist = 1) +
                         netlag(success_list, affirm_adj, pathdist = 1) +
                         #netlag(success_list, correct_adj, pathdist = 2) +
                         netlag(success_list, affirm_adj, pathdist = 2) +
                         #weightlag(success_list, correct_adj) +
                         #   weightlag(success_list, affirm_adj) +
                         # centrality(correct_adj, type = c("indegree")) +
                         centrality(affirm_adj, type = c("indegree")) +
                         covariate(reciprocity_rate_list, coefname = "Reciprocity Rate")+
                         covariate(days_list, coefname = "Days"),
                       #  covariate(age_list) +
                       #  covariate(black_list) +
                       #  covariate(lsiExit_list),
                       re.node = FALSE,
                       time.linear = FALSE,
                       family = binomial(link="logit"))

Affirmations_6 <- tnam(success_list ~
                         #netlag(success_list, correct_adj, pathdist = 1) +
                         netlag(success_list, affirm_adj, pathdist = 1) +
                         #netlag(success_list, correct_adj, pathdist = 2) +
                         netlag(success_list, affirm_adj, pathdist = 2) +
                         #weightlag(success_list, correct_adj) +
                         #   weightlag(success_list, affirm_adj) +
                         # centrality(correct_adj, type = c("indegree")) +
                         centrality(affirm_adj, type = c("indegree")) +
                         covariate(reciprocity_rate_list, coefname = "Reciprocity Rate")+
                         covariate(days_list, coefname = "Days") +
                         covariate(age_list, coefname = "Age"),
                       #  covariate(black_list) +
                       #  covariate(lsiExit_list),
                       re.node = FALSE,
                       time.linear = FALSE,
                       family = binomial(link="logit"))

Affirmations_7 <- tnam(success_list ~
                         #netlag(success_list, correct_adj, pathdist = 1) +
                         netlag(success_list, affirm_adj, pathdist = 1) +
                         #netlag(success_list, correct_adj, pathdist = 2) +
                         netlag(success_list, affirm_adj, pathdist = 2) +
                         #weightlag(success_list, correct_adj) +
                         #   weightlag(success_list, affirm_adj) +
                         # centrality(correct_adj, type = c("indegree")) +
                         centrality(affirm_adj, type = c("indegree")) +
                         covariate(reciprocity_rate_list, coefname = "Reciprocity Rate")+
                         covariate(days_list, coefname = "Days") +
                         covariate(age_list, coefname = "Age") +
                         covariate(black_list, coefname = "Black"),
                       #  covariate(lsiExit_list),
                       re.node = FALSE,
                       time.linear = FALSE,
                       family = binomial(link="logit"))

Affirmations_8 <- tnam(success_list ~
                         #netlag(success_list, correct_adj, pathdist = 1) +
                         netlag(success_list, affirm_adj, pathdist = 1) +
                         #netlag(success_list, correct_adj, pathdist = 2) +
                         netlag(success_list, affirm_adj, pathdist = 2) +
                         #weightlag(success_list, correct_adj) +
                         #   weightlag(success_list, affirm_adj) +
                         # centrality(correct_adj, type = c("indegree")) +
                         centrality(affirm_adj, type = c("indegree")) +
                         covariate(reciprocity_rate_list, coefname = "Reciprocity Rate")+
                         covariate(days_list, coefname = "Days") +
                         covariate(age_list, coefname = "Age") +
                         covariate(black_list, coefname = "Black") +
                         covariate(lsi_list, coefname = "LSI Entry"),
                       re.node = FALSE,
                       time.linear = FALSE,
                       family = binomial(link="logit"))

Affirmations_9 <- tnam(success_list ~
                         #netlag(success_list, correct_adj, pathdist = 1) +
                         # netlag(success_list, affirm_adj, pathdist = 1) +
                         #netlag(success_list, correct_adj, pathdist = 2) +
                         netlag(success_list, affirm_adj, pathdist = 2) +
                         #weightlag(success_list, correct_adj) +
                         weightlag(success_list, affirm_adj) +
                         # centrality(correct_adj, type = c("indegree")) +
                         centrality(affirm_adj, type = c("indegree")) +
                         covariate(reciprocity_rate_list, coefname = "Reciprocity Rate")+
                         covariate(days_list, coefname = "Days") +
                         covariate(age_list, coefname = "Age") +
                         covariate(black_list, coefname = "Black") +
                         covariate(lsi_list, coefname = "LSI Entry"),
                       re.node = FALSE,
                       time.linear = FALSE,
                       family = binomial(link="logit"))

Affirmations_mods <- list(Affirmations_1, Affirmations_2, Affirmations_3, Affirmations_4, Affirmations_5, Affirmations_6, Affirmations_7, Affirmations_8, Affirmations_9)

sink("Table2.txt")
texreg::screenreg(Affirmations_mods,
                  custom.coef.names = c("Intercept",
                                        "Direct Affirmations, Successful Peers",
                                        "Indirect Affirmations, Successful Peers",
                                        "Affirmations Recieved",
                                        "Rate of Affirmations Recieved Reciprocated",
                                        "Days",
                                        "Age",
                                        "Black",
                                        "LSI Entry",
                                        "Weigted Affirmations, Successful Peers"))
sink()


#################
# Interpretation
#################

### Interpretation: Affirmations_8
Affirmations_8noCoef <- tnam(success_list ~
                               #netlag(success_list, correct_adj, pathdist = 1) +
                               netlag(success_list, affirm_adj, pathdist = 1) +
                               #netlag(success_list, correct_adj, pathdist = 2) +
                               netlag(success_list, affirm_adj, pathdist = 2) +
                               #weightlag(success_list, correct_adj) +
                               #   weightlag(success_list, affirm_adj) +
                               # centrality(correct_adj, type = c("indegree")) +
                               centrality(affirm_adj, type = c("indegree")) +
                               covariate(reciprocity_rate_list)+
                               covariate(days_list) +
                               covariate(age_list) +
                               covariate(black_list) +
                               covariate(lsi_list),
                             re.node = FALSE,
                             time.linear = FALSE,
                             family = binomial(link="logit"))

data <- Affirmations_8noCoef$model

newdata.1 <- data.frame(
  netlag.pathdist1 = data$netlag.pathdist1,
  netlag.pathdist2.decay0.5 = mean(data$netlag.pathdist2.decay0.5, na.rm = TRUE),
  indegree = mean(data$indegree, na.rm = TRUE),
  covariate.x = mean(data$covariate.x),
  covariate.y = mean(data$covariate.y),
  covariate.x.1 = mean(data$covariate.x.1),
  covariate.y.1 = median(data$covariate.y.1),
  covariate = mean(data$covariate)
)


predicted.data.1 <- as.data.frame(predict(Affirmations_8noCoef, newdata = newdata.1,
                                          type="link", se=TRUE))

newdata.1 <- cbind(newdata.1, predicted.data.1)

std <- qnorm(0.95 / 2 + 0.5)

newdata.1$ymin <- Affirmations_8$family$linkinv(newdata.1$fit - std * newdata.1$se)
newdata.1$ymax <- Affirmations_8$family$linkinv(newdata.1$fit + std * newdata.1$se)
newdata.1$fit <- Affirmations_8$family$linkinv(newdata.1$fit)

library(ggplot2)
mod1  <- ggplot(data, aes(x=netlag.pathdist1)) +
  geom_ribbon(data = newdata.1, aes(y=fit, ymin=ymin, ymax=ymax), alpha = 0.5) +
  geom_line(data = newdata.1, aes(x = netlag.pathdist1, y=fit), size = 1.5, colour = "firebrick4") +
  scale_y_continuous(limits=c(0.65,1)) +
  theme(legend.position = c(0.2, 0.8),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
  annotate("label", x = 40, y = 0.7, label = 'bold("A: Direct, Successful Peers")', parse = TRUE, size = 6) +
  # ggtitle("A: Direct Affirmations, Successful Peers") +
  labs(x="Successful Alters [Direct Effect]", y="Probability of Graduation")
mod1

interp.data.frame.1 <- newdata.1[1:5,]
interp.data.frame.1$first.order <- NA
interp.data.frame.1[1,]$first.order <- with(newdata.1, mean(netlag.pathdist1) - 2*sd(netlag.pathdist1))
interp.data.frame.1[2,]$first.order <- with(newdata.1,  mean(netlag.pathdist1) - 1*sd(netlag.pathdist1))
interp.data.frame.1[3,]$first.order <- with(newdata.1,  mean(netlag.pathdist1))
interp.data.frame.1[4,]$first.order <- with(newdata.1,  mean(netlag.pathdist1) + 1*sd(netlag.pathdist1))
interp.data.frame.1[5,]$first.order <- with(newdata.1,  mean(netlag.pathdist1) + 2*sd(netlag.pathdist1))

interp.data.frame.1$second.order <- NA
interp.data.frame.1$second.order <-with(newdata.1, mean(netlag.pathdist2.decay0.5))

data$first.order <- data$netlag.pathdist1

first.order.model <- glm(response ~
                           first.order +
                           netlag.pathdist2.decay0.5 +
                           indegree +
                           covariate.x + covariate.y + covariate.x.1 + covariate.y.1 + covariate,
                         data = data,
                         family = binomial)

predicted.data <- as.data.frame(predict(first.order.model, newdata = interp.data.frame.1,
                                        type="link", se=TRUE))

new.data <- cbind(predicted.data)

std <- qnorm(0.95 / 2 + 0.5)
new.data$ymin <- Affirmations_8noCoef$family$linkinv(new.data$fit - std * new.data$se.fit)
new.data$ymax <- Affirmations_8noCoef$family$linkinv(new.data$fit + std * new.data$se.fit)
new.data$fit <- Affirmations_8noCoef$family$linkinv(new.data$fit)

newdata.2 <- data.frame(
  netlag.pathdist1 = mean(data$netlag.pathdist1),
  netlag.pathdist2.decay0.5 = data$netlag.pathdist2.decay0.5,
  indegree = mean(data$indegree, na.rm = TRUE),
  covariate.x = mean(data$covariate.x),
  covariate.y = mean(data$covariate.y),
  covariate.x.1 = mean(data$covariate.x.1),
  covariate.y.1 = median(data$covariate.y.1),
  covariate = mean(data$covariate)
)


predicted.data.2 <- as.data.frame(predict(Affirmations_8noCoef, newdata = newdata.2,
                                          type="link", se=TRUE))

newdata.2 <- cbind(newdata.2, predicted.data.2)
newdata.2$ymin <- Affirmations_8$family$linkinv(newdata.2$fit - std * newdata.2$se)
newdata.2$ymax <- Affirmations_8$family$linkinv(newdata.2$fit + std * newdata.2$se)
newdata.2$fit <- Affirmations_8$family$linkinv(newdata.2$fit)

mod2  <- ggplot(data, aes(x=netlag.pathdist2.decay0.5)) +
  geom_ribbon(data = newdata.2, aes(y=fit, ymin=ymin, ymax=ymax), alpha = 0.5) +
  geom_line(data = newdata.2, aes(x = netlag.pathdist2.decay0.5, y=fit), size = 1.5, colour = "firebrick4") +
  scale_y_continuous(limits=c(0.65,1)) +
  theme(legend.position = c(0.2, 0.8),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))+
  annotate("label", x = 62.5, y = 0.7, label = 'bold("B: Indirect, Successful Peers")', parse = TRUE, size = 6) +
  #ggtitle("B: Indirect Affirmations, Successful Peers") +
  labs(x="Successful Alters [Indirect Effect]", y=NULL)

mod2

interp.data.frame.2 <- newdata.2[1:5,]
interp.data.frame.2$second.order <- NA
interp.data.frame.2[1,]$second.order <- with(newdata.2, mean(netlag.pathdist2.decay0.5) - 2*sd(netlag.pathdist2.decay0.5))
interp.data.frame.2[2,]$second.order <- with(newdata.2,  mean(netlag.pathdist2.decay0.5) - 1*sd(netlag.pathdist2.decay0.5))
interp.data.frame.2[3,]$second.order <- with(newdata.2,  mean(netlag.pathdist2.decay0.5))
interp.data.frame.2[4,]$second.order <- with(newdata.2,  mean(netlag.pathdist2.decay0.5) + 1*sd(netlag.pathdist2.decay0.5))
interp.data.frame.2[5,]$second.order <- with(newdata.2,  mean(netlag.pathdist2.decay0.5) + 2*sd(netlag.pathdist2.decay0.5))

interp.data.frame.2$first.order <- NA
interp.data.frame.2$first.order <-with(newdata.2, mean(netlag.pathdist1))

data$second.order <- data$netlag.pathdist2.decay0.5

second.order.model <- glm(response ~
                            second.order +
                            netlag.pathdist1 +
                            indegree +
                            covariate.x + covariate.y + covariate.x.1 + covariate.y.1 + covariate,
                          data = data,
                          family = binomial)

predicted.data <- as.data.frame(predict(second.order.model, newdata = interp.data.frame.2,
                                        type="link", se=TRUE))

new.data.2 <- cbind(predicted.data)

std <- qnorm(0.95 / 2 + 0.5)
new.data.2$ymin <- Affirmations_8noCoef$family$linkinv(new.data.2$fit - std * new.data.2$se.fit)
new.data.2$ymax <- Affirmations_8noCoef$family$linkinv(new.data.2$fit + std * new.data.2$se.fit)
new.data.2$fit <- Affirmations_8noCoef$family$linkinv(new.data.2$fit)

source("multiplot.R")


png("Figure2.png", width = 900)
multiplot(mod1, mod2, cols = 2)
dev.off()

loocv_predictions_baseline <- list()
loocv_predictions_net1 <- list()
loocv_predictions_net2 <- list()
loocv_predictions_both <- list()

data <- Affirmations_8$model
data <- with(data,
             data.frame(
               success = response,
               net1 = netlag.pathdist1,
               net2 = netlag.pathdist2.decay0.5,
               indegree = indegree,
               reciprocity = `covariate.Reciprocity Rate`,
               days = covariate.Days,
               age = covariate.Age,
               black = covariate.Black,
               lsi = `covariate.LSI Entry`
             ))

colnames(data)

baseline_formula <- as.formula(success ~ indegree + reciprocity + days + age + black + lsi)
net1_formula <- as.formula(success ~ net1 + indegree + reciprocity + days + age + black + lsi)
net2_formula <- as.formula(success ~ net2 + indegree + reciprocity + days + age + black + lsi)
both_formula <- as.formula(success ~ net1 + net2 + indegree + reciprocity + days + age + black + lsi)

for(i in 1:nrow(data)){
  test_data <- data[i,]
  train_data <- data[-i,]
  
  # Base Model
  baseMod <- glm(baseline_formula,
                 model = TRUE,
                 family = binomial(link = "logit"),
                 na.action = na.omit,
                 data = train_data,
                 control = glm.control(epsilon = 0.0001, maxit = 100))
  
  predicted.data.base <- as.data.frame(predict(baseMod, newdata = test_data,
                                               type="response", se=TRUE))
  
  loocv_predictions_baseline[[i]] <- predicted.data.base$fit
  
  # Net 1 Model
  net1Mod <- glm(net1_formula,
                 model = TRUE,
                 family = binomial(link = "logit"),
                 na.action = na.omit,
                 data = train_data,
                 control = glm.control(epsilon = 0.0001, maxit = 100))
  
  predicted.data.net1 <- as.data.frame(predict(net1Mod, newdata = test_data,
                                               type="response", se=TRUE))
  
  loocv_predictions_net1[[i]] <- predicted.data.net1$fit
  
  # Net 2 Model
  net2Mod <- glm(net2_formula,
                 model = TRUE,
                 family = binomial(link = "logit"),
                 na.action = na.omit,
                 data = train_data,
                 control = glm.control(epsilon = 0.0001, maxit = 100))
  
  predicted.data.net2 <- as.data.frame(predict(net2Mod, newdata = test_data,
                                               type="response", se=TRUE))
  
  loocv_predictions_net2[[i]] <- predicted.data.net2$fit
  
  # Both Model
  bothMod <- glm(both_formula,
                 model = TRUE,
                 family = binomial(link = "logit"),
                 na.action = na.omit,
                 data = train_data,
                 control = glm.control(epsilon = 0.0001, maxit = 100))
  
  predicted.data.both <- as.data.frame(predict(bothMod, newdata = test_data,
                                               type="response", se=TRUE))
  
  loocv_predictions_both[[i]] <- predicted.data.both$fit
  
  
  indices <- c(250, 500, 750, 1000, 1250)
  if(i %in% indices){
    print(paste0(i, "iterations complete"))
  }
}

baseMod <- glm(baseline_formula,
               model = TRUE,
               family = binomial(link = "logit"),
               na.action = na.omit,
               data = data,
               control = glm.control(epsilon = 0.0001, maxit = 100))

data$baseMod_fitted <- as.data.frame(predict(baseMod, newdata = data,
                                             type="response", se=TRUE))

net1Mod <- glm(net1_formula,
               model = TRUE,
               family = binomial(link = "logit"),
               na.action = na.omit,
               data = data,
               control = glm.control(epsilon = 0.0001, maxit = 100))

data$net1_fitted <- as.data.frame(predict(net1Mod, newdata = data,
                                          type="response", se=TRUE))

net2Mod <- glm(net2_formula,
               model = TRUE,
               family = binomial(link = "logit"),
               na.action = na.omit,
               data = data,
               control = glm.control(epsilon = 0.0001, maxit = 100))

data$net2_fitted <- as.data.frame(predict(net2Mod, newdata = data,
                                          type="response", se=TRUE))

bothMod <- glm(both_formula,
               model = TRUE,
               family = binomial(link = "logit"),
               na.action = na.omit,
               data = data,
               control = glm.control(epsilon = 0.0001, maxit = 100))

data$both_fitted <- as.data.frame(predict(bothMod, newdata = data,
                                          type="response", se=TRUE))


data$loocv_predictions_baseline <- unlist(loocv_predictions_baseline)
data$loocv_predictions_both <- unlist(loocv_predictions_both)
data$loocv_predictions_net1 <- unlist(loocv_predictions_net1)
data$loocv_predictions_net2 <- unlist(loocv_predictions_net2)

# Comparing PR Curves
library(DMwR)
# v. 0.4.1
library(ROCR)
# v. 1.0-7
library(caTools)
# v. 1.17.1

auc_pr <- function(obs, pred) {
  xx.df <- prediction(pred, obs)
  perf  <- performance(xx.df, "prec", "rec")
  xy    <- data.frame(recall=perf@x.values[[1]], precision=perf@y.values[[1]])
  
  # take out division by 0 for lowest threshold
  xy <- subset(xy, !is.nan(xy$precision))
  
  res   <- trapz(xy$recall, xy$precision)
  res
}

# get plot coordinates
rocdf <- function(pred, obs, data=NULL, type=NULL) {
  # plot_type is "roc" or "pr"
  if (!is.null(data)) {
    pred <- eval(substitute(pred), envir=data)
    obs  <- eval(substitute(obs), envir=data)
  }
  
  rocr_xy <- switch(type, roc=c("tpr", "fpr"), pr=c("prec", "rec"))
  rocr_df <- prediction(pred, obs)
  rocr_pr <- performance(rocr_df, rocr_xy[1], rocr_xy[2])
  xy <- data.frame(rocr_pr@x.values[[1]], rocr_pr@y.values[[1]])
  colnames(xy) <- switch(type, roc=c("tpr", "fpr"), pr=c("rec", "prec"))
  return(xy)
}


# Plots -- full model
# Calculate AUC
bothAUC <- auc_pr(data$success, data$loocv_predictions_both)
# 0.9312698

xyBoth <- rocdf(data$loocv_predictions_both,
                data$success,
                type="pr")

library(ggplot2)
# v 1.0.1
plotBoth <- ggplot(xyBoth, aes(xyBoth$rec, xyBoth$prec)) +
  geom_line(aes(xyBoth$rec, xyBoth$prec), linetype = 1, size = 2, colour = "firebrick4") +
  xlab("Recall") +
  ylab("Precision") +
  annotate("label", x = 0.5, y = 0.683, label = "AUC = 0.9312698") +
  annotate("label", x = 0.5, y = 0.7, label = 'bold("D: Full Model")', parse = TRUE, size = 6) +
  scale_y_continuous(limits=c(0.65,1)) +
  theme(legend.position = c(0.2, 0.8),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))

# Plots -- non-network model
# Calculate AUC
baseAUC <- auc_pr(data$success, data$loocv_predictions_baseline)
# 0.9103609

xyBase <- rocdf(data$loocv_predictions_baseline,
                data$success,
                type="pr")

plotBase <- ggplot(xyBase, aes(xyBase$rec, xyBase$prec)) +
  geom_line(aes(xyBase$rec, xyBase$prec), linetype = 1, size = 2, colour = "firebrick4") +
  xlab("Recall") +
  ylab("Precision") +
  annotate("label", x = 0.5, y = 0.683, label = "AUC = 0.9103609") +
  annotate("label", x = 0.5, y = 0.7, label = 'bold("A: Baseline Model")', parse = TRUE, size = 6) +
  scale_y_continuous(limits=c(0.65,1)) +
  theme(legend.position = c(0.2, 0.8),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))



# Plots -- first-term model
# Calculate AUC
net1AUC <- auc_pr(data$success, data$loocv_predictions_net1)
# 0.9247763

xyNet1 <- rocdf(data$loocv_predictions_net1,
                data$success,
                type="pr")

plotNet1 <- ggplot(xyNet1, aes(xyNet1$rec, xyNet1$prec)) +
  geom_line(aes(xyNet1$rec, xyNet1$prec), linetype = 1, size = 2, colour = "firebrick4") +
  xlab("Recall") +
  ylab("Precision") +
  annotate("label", x = 0.5, y = 0.683, label = "AUC = 0.9247763") +
  annotate("label", x = 0.5, y = 0.7, label = 'bold("B: Direct Effect Model")', parse = TRUE, size = 6) +
  scale_y_continuous(limits=c(0.65,1)) +
  theme(legend.position = c(0.2, 0.8),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))

# Plots -- second-term model
# Calculate AUC
net2AUC <- auc_pr(data$success, data$loocv_predictions_net2)
# 0.9294625

xyNet2 <- rocdf(data$loocv_predictions_net2,
                data$success,
                type="pr")

plotNet2 <- ggplot(xyNet2, aes(xyNet2$rec, xyNet2$prec)) +
  geom_line(aes(xyNet2$rec, xyNet2$prec), linetype = 1, size = 2, colour = "firebrick4") +
  xlab("Recall") +
  ylab("Precision") +
  annotate("label", x = 0.5, y = 0.683, label = "AUC = 0.9294625") +
  annotate("label", x = 0.5, y = 0.7, label = 'bold("C: Indirect Effect Model")', parse = TRUE, size = 6) +
  scale_y_continuous(limits=c(0.65,1)) +
  theme(legend.position = c(0.2, 0.8),
        axis.text=element_text(size=12),
        axis.title=element_text(size=14,face="bold"))
plotNet2

source("multiplot.R")

png("Figure3.png", width = 850, height = 850)
multiplot(plotBase, plotNet2, plotNet1, plotBoth, cols = 2)
dev.off()


#######################
# Description variables
#######################

source("describe.R")


dat <- Affirmations_8$model
dat[dat$covariate.Days < 0,]$covariate.Days <- NA

outcome <- describe(dat$response, type = "binary", plot_title = "TC Graduation")

net1 <- describe(dat$netlag.pathdist1, type = "continuous", plot_title = "Direct Affirmations")

net2 <- describe(dat$netlag.pathdist2.decay0.5, type = "continuous", plot_title = "Indirect Affirmations")

indegree <- describe(dat$indegree, type = "continuous", plot_title = "In-Degree Centrality")

reciprocity <- describe(dat$`covariate.Reciprocity Rate`, type = "continuous", plot_title = "Reciprocity Rate")

days <- describe(dat$covariate.Days, type = "continuous", plot_title = "Days")

age <- describe(dat$covariate.Age, type = "continuous", plot_title = "Age")

black <- describe(dat$covariate.Black, type = "binary", plot_title = "Race")

lsi <- describe(dat$`covariate.LSI Entry`, type = "continuous", plot_title = "Entry LSI")

# rm(list=setdiff(ls(), "dat"))

dat <- Affirmations_9$model
dat[dat$covariate.Days < 0,]$covariate.Days <- NA


weight <- describe(dat$weightlag, type = "continuous", plot_title = "Weighted Affirmations")

mean(dat$weightlag)
sd(dat$weightlag)


dist <- list(outcome = outcome$dist, direct = net1$dist, indirect = net2$dist,
             indegree = indegree$dist, reciprocity = reciprocity$dist, days= days$dist,
             age = age$dist, black = black$dist, lsi = lsi$dist, weight = weight$dist)

sink("Table1.txt")
dist
sink()

png("Figure1.png", width = 700, height = 700)
multiplot(outcome$summary_plot, net1$summary_plot, net2$summary_plot,
          indegree$summary_plot, reciprocity$summary_plot, days$summary_plot,
          age$summary_plot, black$summary_plot, lsi$summary_plot, weight$summary_plot, cols=2)
dev.off()